!
! astrometry
!
!
! Copyright © 2011-2 F.Hroch (hroch@physics.muni.cz)
!
! This file is part of Munipack.
!
! Munipack is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
! 
! Munipack is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
! 
! You should have received a copy of the GNU General Public License
! along with Munipack.  If not, see <http://www.gnu.org/licenses/>.
!

! TODO:
!
!   * add all WCS related projections
!   * add catalogue/reference matching identifier
!
!  testing:
!   * improve high-precision algorithm to fit all parameters simulatenously
!   * add least squares algorithm
!   * add tunable parameter for units in FITS headers (uncertainitilies
!     + S0): ",',degrees


program astrometer

  use astrofits
  use astrofitting
  use astromatcher

  implicit none

  integer, parameter :: dbl = selected_real_kind(15)
  real(dbl), parameter :: rad = 57.295779513082322865_dbl
  integer, parameter :: clen = 80
  character(len=32768) :: line, record, ref = '', cat = '', output = '', filename = ''
  character(len=clen) :: key_alpha = 'RA', key_delta = 'DE'

  character(len=clen) :: type = 'GNOMONIC', method = 'ROBUST', aunits = '' 
  character(len=clen) :: rframe = '', idcat = ''
  real(dbl), allocatable, dimension(:) :: alpha, delta, x, y, ares, dres
  integer,  allocatable, dimension(:) :: id1,id2
  integer,  allocatable, dimension(:) :: i1,i2
  character(len=clen),allocatable, dimension(:) :: com
  real(dbl), dimension(2) :: crpix, crval, dcrval
  real(dbl) :: s0,sc,dsc,pa,dpa,rms,like,a,da,b,db,s,ds,w,h,wref,href, sm
  real(dbl) :: validity = 0.5, tol = 10.0
  integer :: cc,eq,n,ndat,istat
  integer :: nmatch = 5, maxmatch = 30
  logical :: fcrpix = .false., fcrval = .false., &
       fscale = .false., fangle = .false., matching = .true., initpar = .true.,&
       verbose = .false., wcssave = .false.


  ndat = 0
  n = 0

  ! read input parameters
  do
     read(*,'(a)',end=99,iostat=istat) record
     if( istat /= 0 ) then
        write(*,*) 'Unexpected string on input.'
        stop 21
     end if

     write(*,*) trim(record)


     cc = index(record,'#')
     if( cc == 0 ) then
        line = record
     else
        line = record(1:cc-1)
     end if

     eq = index(line,'=')

    ! astrometrical parameters

     if( index(line,'PROJECTION') > 0 .and. eq > 0) then

        ! IMPORTANT: the projection must * PRECEDE * both CAT and REF

        read(line(eq+1:),*,iostat=istat) type
        if( istat /= 0 ) then
           write(*,*) 'PROJECTION string read failed.'
           stop 21
        end if
        
     else if( index(line,'METHOD') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) method
        if( istat /= 0 ) then
           write(*,*) 'METHOD string read failed.'
           stop 21
        end if

     else if( index(line,'AUNITS') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) aunits
        if( istat /= 0 ) then
           write(*,*) 'AUNITS string read failed.'
           stop 21
        end if

     else if( index(line,'VERBOSE') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) verbose
        if( istat /= 0 ) then
           write(*,*) 'VERBOSE parameter read failed.'
           stop 21
        end if

     else if( index(line,'WCSSAVE') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) wcssave
        if( istat /= 0 ) then
           write(*,*) 'WCSSAVE parameter read failed.'
           stop 21
        end if

     else if( index(line,'CRPIX') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) crpix
        if( istat /= 0 ) then
           write(*,*) 'CRPIX parameter read failed.'
           stop 21
        end if
        fcrpix = .true.

     else if( index(line,'CRVAL') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) crval
        if( istat /= 0 ) then
           write(*,*) 'CRVAL parameter read failed.'
           stop 21
        end if
        fcrval = .true.

     else if( index(line,'SCALE') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) sc
        if( istat /= 0 ) then
           write(*,*) 'SCALE parameter read failed.'
           stop 21
        end if
        fscale = .true.

     else if( index(line,'ANGLE') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) pa
        if( istat /= 0 ) then
           write(*,*) 'ANGLE parameter read failed.'
           stop 21
        end if
        fangle = .true.

     else if( index(line,'TOL') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) tol
        if( istat /= 0 ) then
           write(*,*) 'TOL parameter read failed.'
           stop 21
        end if
        fangle = .true.

     else if( index(line,'INITPAR') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) initpar
        if( istat /= 0 ) then
           write(*,*) 'INITPAR parameter read failed.'
           stop 21
        end if


     ! match parameters

     else if( index(line,'MATCHING') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) matching
        if( istat /= 0 ) then
           write(*,*) 'MATCHING parameter read failed.'
           stop 21
        end if

     else if( index(line,'NMATCH') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) nmatch
        if( istat /= 0 ) then
           write(*,*) 'NMATCH parameter read failed.'
           stop 21
        end if

     else if( index(line,'MAXMATCH') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) maxmatch
        if( istat /= 0 ) then
           write(*,*) 'MAXMATCH parameter read failed.'
           stop 21
        end if

     else if( index(line,'VALIDITY') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) validity
        if( istat /= 0 ) then
           write(*,*) 'VALIDITY parameter read failed.'
           stop 21
        end if


     ! files to working on

!     else if( index(line,'OUTPUT') > 0 .and. eq > 0) then
!
!        read(line(eq+1:),*,iostat=istat) output
!        if( istat /= 0 ) then
!           write(*,*) 'OUTPUT parameter read failed.'
!           stop 21
!        end if

     else if( index(line,'KEY_ALPHA') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) key_alpha
        if( istat /= 0 ) then
           write(*,*) 'KEY ALPHA parameter read failed.'
           stop 21
        end if

     else if( index(line,'KEY_DELTA') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) key_delta
        if( istat /= 0 ) then
           write(*,*) 'KEY DELTA parameter read failed.'
           stop 21
        end if

     else if( index(line,'REF') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) ref
        if( istat /= 0 ) then
           write(*,*) 'REF parameter read failed.'
           stop 21
        end if

        call readfile(ref,wref,href,alpha,delta)
        if( initpar ) call init(type,alpha,delta,crval(1),crval(2),sc)

     else if( index(line,'CAT') > 0 .and. eq > 0) then

        read(line(eq+1:),*,iostat=istat) cat
        if( istat /= 0 ) then
           write(*,*) 'CAT parameter read failed.'
           stop 21
        end if

        call readtable(cat,key_alpha,key_delta,alpha,delta,rframe,idcat)
        if( initpar ) call init(type,alpha,delta,crval(1),crval(2),sc)

     else if( index(line,'FILE') > 0 .and. eq > 0) then

        if( cat == '' .and. ref == '' ) then
           write(*,*) 'Catalogue or reference image unspecified.'
           stop 21
        end if

        read(line(eq+1:),*,iostat=istat) filename, output
        if( istat /=0 ) then
           write(*,*) 'Invalid filename string.'
           stop 21
        end if

        write(*,'(a)') "=> Astrometry calibration for `" // trim(filename) // "'."
        write(*,'(a)') "=> Processing..."

        if( fcrpix ) then
           call readfile(filename,w,h,x,y)
        else
           call readfile(filename,w,h,x,y,crpix)
        end if

        if( initpar ) sc = sc/(sum(crpix)/size(crpix))
!        write(*,*) crval1,crval2,sc

        write(*,'(a)') "=> Matching ..."
        if( matching ) then

           write(*,*) alpha(1),delta(1),x(1),y(1)

           call astromatch(nmatch,maxmatch,type,alpha,delta,x,y, &
                crval(1),crval(2),crpix(1),crpix(2),w,h,id1,id2,like)

        else

           ! tol = tol/3600 in degs????
           call listmatch(type,alpha,delta,x,y,crval(1),crval(2),crpix(1),crpix(2), &
                w,h,sc,pa,tol,id1,id2)
           like = 2.0*validity

        end if

        if( like < validity ) then

           write(*,*) "Mutual match for files `",trim(filename),"' and `", &
                trim(cat),trim(ref),"' failed. Validity under limit: like=", &
                like,' < validity=',validity           
        else

           write(*,'(a)') "=> Fitting ..."
           call astrometry(method,type,alpha,delta,x,y,w,h,crpix(1),crpix(2), &
                crval(1),dcrval(1),crval(2),dcrval(2),sc,dsc,s0,rms,pa,dpa,ares,dres,&
                id1,id2,i1,i2)

           ! use a return flag instead ?
           if( size(i1) > 0 ) then

              allocate(com(13+size(i1)))
              
              a = cos(pa/rad)
              da = (dpa/rad)*sin(pa/rad)
              b = sin(pa/rad)
              db = (dpa/rad)*cos(pa/rad)
              s = 3600.0_dbl*sc
              ds = s*dsc/sc

              ! add catalogue ID !!!!!
              
              if( aunits == '' ) then
                 if( log10(3600.0*rms) < 2.0 ) then
                    aunits = 'arcsec'
                 else if( log10(60.0*rms) < 2.0 ) then
                    aunits = 'arcmin'
                 else 
                    aunits = 'deg'
                 end if
              end if

              write(com(1),'(a)') "Astrometric Solution:"
              write(com(2),'(a,i0)') "Number of objects used = ",size(i1)
              if( aunits == 'arcsec' ) then
                 sm = 3600.0_dbl
                 s = sm/sc
                 ds = s*dsc/sc
                 write(com(3),'(a,g20.10,a)') "Residual sum = ",sm**2*s0,' ["]'
                 write(com(4),'(a,g20.10,a)') "RMS = ",sm*rms,' ["]'
                 write(com(5),'(a,g20.10,a,es9.2,a)') "Scale = ",s," +- ",ds,' ["/pix]'
              else if( aunits == "arcmin" ) then
                 sm = 60.0_dbl
                 s = sm/sc
                 ds = s*dsc/sc
                 write(com(3),'(a,g20.10,a)') "Residual sum = ",sm**2*s0," [']"
                 write(com(4),'(a,g20.10,a)') "RMS = ",sm*rms," [']"
                 write(com(5),'(a,g20.10,a,es9.2,a)') "Scale = ",s," +- ",ds," ['/pix]"
              else if( aunits == "deg" ) then
                 sm = 1.0_dbl
                 write(com(3),'(a,g20.10,a)') "Residual sum = ",s0,' [deg]'
                 write(com(4),'(a,g20.10,a)') "RMS = ",rms,' [deg]'
                 write(com(5),'(a,g20.10,a,es9.2,a)') "Scale = ",sc," +- ",dsc,' [deg/pix]'
              end if
              write(com(6),'(a,f15.10,a,es9.2)') "cos(pa) = ",a," +- ",da
              write(com(7),'(a,f15.10,a,es9.2)') "sin(pa) = ",b," +- ",db
              write(com(8),'(a,f15.10,a,es9.2,a)') "Position Angle (pa)  = ",pa," +- ",dpa," [deg]"
              write(com(9),'(a,f15.10,a,es9.2,a)') "Alpha center projection (CRVAL1) = ", &
                   crval(1)," +- ",dcrval(1)," [deg]"
              write(com(10),'(a,f15.10,a,es9.2,a)') "Delta center projection (CRVAL2) = ", &
                   crval(2)," +- ",dcrval(2)," [deg]"
              write(com(11),'(a,f9.3,a)') "Horizontal center (CRPIX1) = ",crpix(1),' [pix]'
              write(com(12),'(a,f9.3,a)') "Vertical   center (CRPIX2) = ",crpix(2),' [pix]'

              if( aunits == 'arcsec' ) then
                 write(com(13),'(a)') 'Catalogue RA,DEC [deg]      Objects X,Y [pix]     Residuals["]'
              else if( aunits == 'arcmin' ) then
                 write(com(13),'(a)') "Catalogue RA,DEC [deg]      Objects X,Y [pix]     Residuals[']"
              else if( aunits == 'deg' ) then
                 write(com(13),'(a)') 'Catalogue RA,DEC [deg]      Objects X,Y [pix]     Residuals[deg]'
              end if
              

              do n = 1, size(ares)
                 write(com(13+n),'(2f13.8,2f9.3,2en12.1)') &
                      alpha(i1(n)),delta(i1(n)),x(i2(n)),y(i2(n)),sm*ares(n),sm*dres(n)
              end do

              write(*,'("=C> ",a)') (com(n),n=1,size(com))

              
              if( wcssave ) then
                 write(*,'(a)') "=> Saving ..."
                 call wcsupdate(filename,output,type,crval,dcrval,sc,dsc,pa,dpa,fcrpix,crpix,com)
              end if

              deallocate(com)
              write(*,'(a)') "=> Successful."

           else
              write(*,'(a)') "=> Astrometry of `" // trim(filename) // "' unsuccessful."
              write(*,'(a)') "=> Failed."
           end if
           
           deallocate(ares,dres,i1,i2)

        end if

        deallocate(x,y)
        deallocate(id1,id2)

     end if

  end do
99 continue

  if( allocated(alpha) ) deallocate(alpha,delta)

  stop 0


end program astrometer
